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W~) ■ Abstract 

The inability of standard non-interacting cold dark matter (CDM) to account for the 
small scale structure of individual galaxies has led to the suggestion that the dark 
matter may undergo elastic and/or inelastic scattering. We simulate the evolution 
Qh' of an isolated dark matter halo which undergoes both scattering and annihilation. 

Annihilations produce a core that grows with time due to adiabatic expansion of 
the core as the relativistic annihilation products flow out of the core, lessening 
the binding energy. An effective annihilation cross section per unit mass equal to 
.03 cm 2 g _1 (100 km s _1 /v) with a scattering cross section per unit mass of .6 cm g _1 
produces a 3 kpc core in a 10 10 M halo that persists for 100 dynamical times. The 
same cross section leads to a core of only 120 pc in a rich cluster. In addition to 
creating to cores, annihilation should erase structure on scales below ~ 3x 10 8 M . 
Annihilating dark matter provides a mechanism for solving some of the problems 
of non-interacting CDM, at the expense of introducing a contrived particle physics 
model. 
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1 Introduction 



Discrepancies between the Cold Dark Matter (CDM) model and observations 



on galactic scales ( 


Salucci, 2001 




Vloore, 1994|; [Flores & Primack, 1994; 


Burk- 


ert, 1995; 


Salucci & Burkert, 200( 


); Moore, 2001; 


Sellwood & Kosowsky, 2000]; 




Weiner et al., 2001b) have led to recent revival of the suggestion that the dark 



matter is self-interacting flSpergel fc Steinhardt, 2000]) . 



Simulations of halo formation in a CDM universe predict that dark matter 
halos have a cuspy inner profile. Most simulations find that, at small radii, 
p oc r a , with a ~ — 1 to —1.5 ( [Navarro et al., 1996b| ; |Chigna et al., 2000| ; 
Jing fc Suto, 



2000 



Klypin et al., 2000| ). An inner density profile this steep 
may be consistent with the dark matter density profile inferred from many 
high surface brightness galaxies and clusters of galaxies ([Eke et al. (2000| ). 
but see [Borriello fc Salucci (2001| ) ). However, low surface brightness galaxies, 
and dwarf irregular galaxies in particular, require an inner density profile 
which is less steep flMoore (TgggP ; |Burkert (19"9"5| ); [b'lores & Primack (Tggg) ; 
for a contradictory view see van den Bosch fc Swaters (2000] )). In addition, 
Debattista fc Sellwood (2000|) have emphasized that too much dark matter in 
the centers of galaxies is incompatible with rotating bars, which are common 
in massive galaxies such as the Milky Way. 



A separate problem for standard CDM model is that it predicts many more 
dark matter satellites of a Milky Way-sized galaxy than are observed ( [Moore 
et al., 1999| ; [Klypin et al., 1999| ). This substructure cannot simply be dark 
satellites in which stars have not formed, because they would thicken the disk 
of the Milky Way as they passed through the plane of the disk. 



While it is possible that the inclusion of baryonic physics would lead to a 
less cuspy dark matter halo, perhaps due to supernova driven gas outflow in 
dwarf galaxies ( [Navarro et al., 1996a| ), this seems unlikely, as cooling baryons 
would lead to further concentration of the dark matter halo through adiabatic 
compression. Exploding gaseous disks will drive a strong wind off the plane, 
but most simulations show that the cold, dense component of a multi-phase 
medium would be difficult to remove ( [Mac how, 200 1| ). Furthermore, these 
solutions to the cusp problem does not address the problem of too much 
substructure on small scales. 



Recent attention has focused instead on the possibility that the dark matter is 
self-interacting. A remarkable number of different models have been proposed. 
For example, several authors have considered a model in which the dark matter 
is a self-interacting scalar field with quartic potential ( [Peebles, 2000| ; [Peeblesj 
fc Vilenkin, T999| ; |Coodman, 20001 ; |Riotto fc Tkachev, 200q) . [Hu et al. (20001) 



have suggested a model in which cores are produced in dwarf galaxies due to 
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the wave nature of dark matter, if the matter is an ultra-light scalar particle. 
Another interacting dark matter candidate is Q-balls, non-topological soli- 
tons with properties that allow it to scatter, possibly with a time-dependent 
cross-section, or to merge ([Kusenko fc Steinhardt, 20011) . Several authors have 
simulated halo evolution and formation in scenarios in which the dark matter 
scatters. All authors agree that on short time scales, scattering reduces the 
logarithmic slope of the inner density profile substantially; some authors have 
found that the inner density profile subsequently steepens as the halo un- 
dergoes core collapse ( [Kochanek fc White, 2U(J(J| ; |Yoshida et al., 20TIUD , while 



others disagree ( [Burkert, 2000| ; pave et al., 2001| ) 



Self-annihilating dark matter is another possible solution to the problem QKapling- 
hat et al., 2000| ), although the particle physics must be carefully contrived so 
that the dark matter does not completely self-annihilate in the early universe. 
An analytic calculation of the density of a self-annihilating dark matter halo 
demonstrates the formation of a core. Physically, the core develops because 
annihilation becomes very efficient at removing mass from the halo above 
a density that depends on the annihilation cross section. |Kaplinghat et al. 



(2000|) point out that this annihilation also causes an adiabatic expansion of 
the halo, and numerically calculate the effects of adiabatic expansion; the core 
density including expansion is smaller than that predicted from annihilation 
alone by a factor of ten (for dwarf galaxies) to three (for clusters). In addi- 
tion, [Hui (2001"]) points out that unitarity bounds require that annihilation be 
accompanied by scattering. 

In this paper, we simulate the behavior of an isolated dark matter halo which 
undergoes both scattering and annihilation. In section 2, we describe our dark 
matter model and our simulation methods. In section 3 we present the results 
of our simulations, and in section 4 we compare our results with those of 
previous authors. 



2 Model 



2. 1 Initial halo and units 



We simulated the evolution of an isolated halo; the initial density profile is 
the [Hcrnquist (1990| ) profile, 



M 

2nr(r + r s 



P= n--,_ , - N3 » W 
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where M is the total mass of the halo, and r s is the scale radius. We take the 
dynamical time for this model to be the circular orbit time at the scale radius 



yn 



4ir\ 



GM 



(2) 



In the simulations below, we work in units in which G — 1, M — 1, and r s — 1. 

The profile has been generated to guarantee that the center of mass of the halo 
is also the origin used in calculating the density of the Hernquist profile. We 
first choose the energy for each particle, then select a radius for the particle; 
these determine the magnitude of the velocity of the particle. We adopt an 
isotropic velocity distribution. For each particle generated with position r, 
a second particle is generated with position — r; the second particle has the 
same energy as the first. The distribution is truncated at a radius of IOOtv 
The direction of the velocity of each particle is generated independently, but 
the net velocity of the halo is set to zero. 



2.2 Model for scattering and annihilation. 



We adopt a model similar to that used by |Burkert (2000| ). The probability for 
interaction for a particle which moves a distance Ax in time step At is 



where A = l/(no~) is the mean free path for interactions; this expression is valid 
only when Ax <C A. Our simulations include both scattering and annihilation. 
We use the notation a* to denote the annihilation cross section per unit mass 
and a* to denote the scattering cross section per unit mass. 

The probability that a pair of particles scatters in time step At is 

P s ,ij = o-* s VijPiAt, (4) 



where p\ is the local density, and is the relative velocity of particles i and 
j. Scattering is completely elastic, and is isotropic in the center of momentum 
frame of the pair of particles. 

As pointed out by [Kaplinghat et al. (2000| ) , only s-wave annihilation is consis- 
tent with observations. For s-wave annihilation, o~\v\ is independent of velocity, 
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so that the cross section for annihilation per unit mass, <r*, is given by 

* °~a V a , v 

°a = , ( 5 ) 

m v 



where a a v a is determined by the scattering potential. The numerical value we 
assume for the scattering cross section is given below. 

The probability for an annihilation interaction in time step At is: 

P a = a* aPl vAt = ^ Pl At, (6) 



where pi is the local density. The annihilation probability is independent of 
velocity. 

We model the decay of the dark matter particles represented by our simulation 
particles by assigning to each simulation particle a probability for annihilation 
given by equation (6). If the local density is large enough, P a is the probability 
that some of the dark matter particles represented by the simulation particle 
would interact and decay. If a simulation particle decays, its mass rrii is reduced 
by fraction F, 

m i:new = (1 - F)m ij0ld . (7) 



Thus, the effective probability for annihilation of individual dark matter par- 
ticles in our model is the product FP a . The mass lost in an annihilation is 
a pre-set fraction of the current mass of the particle undergoing annihilation. 
In other words, the fractional mass lost in an annihilation is held fixed. The 
absolute mass loss, m ij0 i d — m ijnew , decreases as a particle undergoing repeated 
annihilation loses mass. 



2.3 Numerical Parameters 



Our simulations were run on a GRAPE3-AF special purpose hardware card 
QOkumura et al., 1993|) using a direct summation code. The time step was cho- 
sen to be small enough to satisfy both the Courant condition and to ensure 
that vAt < A. Each run contained 30,000 particles. Detailed numerical pa- 
rameters for each run are shown in Table 1. Runs S0A0, S3A0, and S3A1 ran 
for roughly 10 dynamical times with approximately 2100 steps per dynamical 
time; a fourth run, S.3A.1, was done with a much larger Plummer softening 
and time step; it ran for 100 dynamical times with approximately 400 steps 
per dynamical time. 
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Table 1 

Simulation Parameters 



Name 






e 


At 


A steps 


if/tdyn 


SOAO 








0.01 


0.006 


20,000 


9.55 


S3A0 


3.0 





0.01 


0.006 


15,750 


7.52 


S3A1 


3.0 


1.5 


0.01 


0.006 


20,000 


9.55 


S.3A.1 


0.3 


0.15 


0.05 


0.03 


40000 


95.5 



One run, SOAO, had no interactions, and was used to determine the stability 
of the integration. The run S3A0 had no annihilations, but non-zero scattering 
cross-section. Two runs with annihilation and scattering were performed. One, 
S3A1, was short (the final time was a tenth of a Hubble time assuming a 
1O 1O M halo), but had the same scattering cross-section as the scattering-only 
run. The other, S.3A.1, ran for ten times longer, with both annihilation and 
scattering cross section reduced by a factor of ten. In both annihilation runs 
the scattering and annihilation probabilities are equal at the scale radius; inside 
of that radius scattering will be enhanced relative to annihilation because the 
scattering probability increases with velocity. The fractional mass loss per 
annihilation in both was F = 5%. 

The cross section per unit mass for scattering can be calculated in physical 
units from 

^ = A ^[-^f-){^- c ) cmg . (8) 



The cross-section for annihilation is given by 

as * f 1Ql ° M rA ( R s V /I00kms^\ 2 _ x 



For R s = 2 kpc and M = 10 10 M Q , the half mass dynamical time is .168 Gyr. 
Our cross-sections are a* = 5.76cm 2 g _1 and a* = 5.88cm 2 g _1 (100kms _1 /t') 
for the run S3A1; in the run S.3A.1, both cross sections are a factor of ten 
smaller. Our values of a* span the range suggested by [Wandeft et al. (2000| ) . 



Kochanek fc White (2000]) considered a broader range of values for a*; our 



choice for our runs with the largest scattering cross section is in the middle of 
their range. 

In both annihilating runs, we chose our annihilation cross section to be com- 
parable to the scattering cross section. However, this annihilation cross section 
in the run S3A1 is about 300 times larger than that suggested by [Kapling 



hat et al. (2000[) based on fits of rotation curves to dwarf galaxies. The run 
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S.3A.1, with an annihilation cross section ten times smaller, is closer to what 
they suggest. With a mass loss factor F of 5%, the effective cross section Fa* 
is comparable to the value required to match observations. In both cases, the 
scattering cross section is above the minimum required by unitarity bounds 
QHui, 2001|). 



3 Results 



The inclusion of annihilation in our simulations causes a reduction in the 
central density of the halo. Figure 1 shows the evolution in the density profile 
of a halo with both annihilation and scattering; the top panel shows a short 
run with high cross section (S3A1), and the bottom a run with smaller cross 
section run for longer time (S.3A.1). The central density rapidly decreases, 
and remains low for the duration of the simulations, about 100 (10) dynamical 
times for run S.3A.1 (S3A1). 

The core is not due to numerical effects. The force softening is .05 (.01) in 
the run S.3A.1 (S3A1); the flattening of the inner density profile we observe 
occurs at radii 2-10 (10-50) times larger than that and is unlikely to have been 
caused by softening. To test the stability of our code and the effects of our 
relatively limited particle number, we ran a simulation with no interactions 
for ten dynamical times. The simulation with no interactions S0A0 shows no 
evidence of the development of a core. The density profile of the halo was 
essentially unchanged over the course of the simulation. Figure 1 shows the 
beginning profile and the final profile for this run. 

We have examined the mechanism by which the core forms in our interaction 
simulations. There are three related mechanisms one might imagine contribut- 
ing to the development of a core. First, mass lost due to annihilation and the 
disappearance of the relativistic decay products from the halo directly de- 
creases the central density of the halo. Second, this direct mass loss reduces 
the gravitational potential of the halo, causing particles in the central region 
to become less bound. Third, energy is transferred from the central region of 
the halo outward by scattering, depressing the central density. 

Although scattering at first suppresses the central density, core collapse begins 
fairly rapidly in halos with scattering only, as reported by |Kochanek fc White 



(2000 ). We have run one simulation with scattering alone to determine whether 
the core that initially forms persists. The results of this run are shown in 
Figure 2. The central density initially drops and then increases after a few 
dynamical times.Q In the absence of annihilation, scattering does not produce 



This run was terminated after 5.8 dynamical times. At that point, the central 
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10.00 F 




Fig. 1. Evolution of density profile in two models with annihilation and scattering. 
The heavy lines in both plots show the chgnge in core behavior without interactions. 
In the top plot, the force softening length is .01; in the bottom plot it is 0.05. 
In this and later density plots, only the central region of the halo is shown; the 



a persistent core in simulations of isolated halos. It is unclear whether or not 
cores persist in high resolution simulations of halo formation in a cosmological 
setting QYoshida et al., 2000| ; pave et al., 2001] ) . 




Fig. 2. Evolution of density profile in the scattering-only model S3A0. The core 
density initially decreases, but then begins to grow with time. 

In order to differentiate between direct mass loss from annihilation and a 
change in the density profile due to a change in gravitational potential, we 
calculate the change in number density of simulation particles as a function 
of radius. If direct mass loss due to annihilation is causing most of the decline 
seen in the density profile, the number density profile ought to be unchanged, 



density became so large that the typical particle was moving more than one mean 
free path in a time step. 
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with particles in the central region of the halo simply becoming less massive. 
Instead, as can be seen in Figure 3, which shows the evolution in the number 
density of particles early in the simulation, we find that the number density 
rapidly diminishes in the center of the halo, indicating that much of the change 
in central density is due to the motion of particles away from the center, not 
from direct mass loss. The trend of particle number density decreasing in the 
center continues throughout the simulation, as can be seen from Figure 1. 
Density profiles are calculated with a fixed number of particles per bin, so 
the location of the first point in the density profile is a radius that encloses a 
fixed number of particles; this point moves farther from the origin throughout 
the simulation. Furthermore, the total mass lost by the halo over the entire 
simulation is only 12% of the initial halo mass. 




Fig. 3. Evolution in number density of particles in annihilation+scattering model. 
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The halo in the annihilating model forms a core very rapidly, and the size of the 
core continues to grow over the course of the simulation. We fit two different 
functional forms to the density profile of the halo. One is a modification of 
the initial Hernquist profile 1 to allow for a core, 

PC {x + x^l + x) 3 ' 1 ' 



where x = r/r s , and x c = r c /r s is the core radius in units of the scale radius. 
The other density profile was derived by [Kaplinghat et al. (2000| ) ; their profile 
takes into account modifications of the initial halo density profile due to direct 
annihilation only, without scattering. Adiabatic expansion of the core as mass 
is lost decreases the core density. We have derived the profile appropriate for 
an initial Hernquist model, 

_ P0 x c / -i -i \ 

P ~ x(l + xf + x c 1 ' 



(note that x c is related to the annihilation cross section in Kaplinghat 's model 
(Eq. 13); in doing our fits we have treated free parameter to see how well 

the shape of the density profile is fit by the function). Fits to both functions 
are shown at four different epochs in Figure 4. The two profiles differ only in 
the innermost region of the halo. Burkert (1995|) found that dwarf galaxies are 
well fit by the profile 

P B ~ 7 1 \( 2 i 2V ^ > 

(r + r c )[r z + r z ) 



where po is the central density of the halo, and r c its core radius. We do not 
expect the Burkert profile to be a good fit to our halo because the initial halo 
profile has both a different radial dependence at large radii, and a length scale 
(the scale radius) in the initial conditions separate from the core radius. 

The halo shows a clear evolution from a Kaplinghat profile to a modified 
Hernquist profile. At early times, the Kaplinghat model is a better fit than the 
modified Hernquist model. At later times, the modified Hernquist model is a 
better fit. This is good news, because the modified Hernquist profile is closer 
in form to the Burkert profile. The evolution in profile shape makes sense; 
early in the simulation the adiabatic effect of annihilation on the gravitational 
potential will be small, and the Kaplinghat profile should be a good fit. The 
core radius in the modified Hernquist profile continues to grow throughout the 
simulation, as shown in Figure 5. The final core radius is 3 kpc, in the range 
suggested by Purkert (1995Q for a halo of the mass we studied. 
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0.1 



1.0 

r/r c 



0.1 



1.0 

r/r c 



10.0 



Fig. 4. Fits to the density profile of the annihilating dark matter run S.3A.1. The 
solid line shows the best fitting modified Hernquist profile, the dotted line the 
best-fitting Kaplinghat profile, and the data are the symbols. As the simulation 
proceeds, the halo is better fit by a modified Hernquist profile than a Kaplinghat 
profile. 



4 Discussion 



Our simulations demonstrate that annihilation produce cores of radius a few 
kiloparsecs in dark matter halos of dwarf galaxies; we now discuss whether 
this dark matter model preserves the successes of non- interacting CDM, the 
extent to which it addresses the problem of too much substructure in galactic 
dark matter halos, and the plausibility of the particle physics of the model. 
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Fig. 5. Growth of the core radius with time in the annihilating DM model 

There is little evidence for cores in clusters; for an annihilation model to be 
plausible it must modify small scale behavior of dark matter while leaving 
large scale behavior unchanged. In the analytic model of [Kaplinghat et al. 
(2000|) , the core radius grows linearly in time, and is given by 

T 

x c = — = p s a a t/m = p s Fa* a t (13) 



where p s is the characteristic density of the NFW profile and t is the age of the 
halo. Since lower mass halos have higher density, and formed earlier, than high 
mass halos, one would expect a decreasing core radius with increasing mass. 
Figure 6 shows the trend of core radius with mass, assuming the core radius 
is given by equation (13); the core radius has been normalized to that of a 
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10 10 h~ 1 M Q halo (approximately 3 kpc) to eliminate the dependence of the core 
radius on cross section. Characteristic densities and halo ages were calculated 
using the [Eke et al. (2000[ ) model, assuming a Q = .3, = 1 — = .7 
flat cosmology with h = .65. Most of the decrease in core radius is due to 
the decrease in characteristic density as mass increases; the age difference of 
the lowest and highest mass halo accounts for only 20% of the difference in 
core radius. Our model clearly creates large cores in low mass/high density 
halos without producing them in high mass/low density halos. Thus if the 
annihilations lead to a core radius of 3 kpc for dwarf galaxies, as appropriate, 
the predicted dark matter core radius in rich clusters would be ~ 120 pc, 
which is deep within the central galaxy. 



1.000 



0.100 
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10'° 10' 
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Fig. 6. Relative size of cores radius in halos of different mass. 
The core structure of the dark matter halo of a high surface brightness galaxies 
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is unclear. |Eke et al. (2000|) argue that the concentrations of dark matter halos 



predicted in the standard ACDM model are consistent with both the Tully- 
Fisher relationship and the amount of dark matter interior to the solar radius. 
On the other hand, [Weiner et al. (200ia| ) argues that detailed mass modeling 
of the barred spiral NGC 4123 requires either a modest core in the dark 
matter halo, with a core radius of 1-6 kpc, or an NFW profile with very low 
concentration. Arguing more generally, pebattista fc Sellwood (200(J| , |1998| ) 
find that a high density dark matter halo in a barred galaxy leads to rapid 
slowing of bar rotation. For the same choice of annihilation cross-section that 
generates a 3 kpc core in dwarf galaxies, the core radius in a 10 12 M galaxy 
such as the Milky Way would be ~ 1.2 kpc, and would resolve the difficulties 
of bar spirals in standard ACDM discussed by |Debattista fc Sellwood (2000| ). 



The standard ACDM model also predicts more substructure on galactic scales 
than is observed ( |Moore et al., 1999] ; [Klypin et al., 199"9|) . Although the sim- 



ulations we have done cannot directly address the question of substructure, 
Figure 6 can be extended to low mass halos. As halo mass shrinks, the core 
radius increases while the virial radius increases. The two are equal for halos 
of mass 3xlO 7 M ; halos with mass 10 8 solar masses have a core radius that is 
half the virial radius. This suggests that halos smaller than a few times 10 s so- 
lar masses would not be present in an annihilating dark matter model. Further 
simulation of structure formation in a cosmological context with annihilating 
dark matter would be helpful. 

One of the problems with earlier self-interacting dark matter proposals ( |Carl-| 
son et al., 1992| Machacek, 1994|) is that ram pressure stripping will destroy 



small halos too efficiently ( pe Laix et al., 1995| ). Newer models in which 
the dark matter behaves like a viscosity-free fluid ( Peebles, 2000| ; Peebles fc" 



Vilcnkim 1999j ; |Goodman, 200C ; Riotto fc Tkachev, 2000 ) avoid this problem 



but may overly- suppress the growth of structure on small scales in the linear 
regime ( Peebles, 2000|) , a problem shared by the earlier self-interacting dark 



matter models ( |Machacek, 1994| ). 



The model proposed by |Carlson et al. (1992| ) has some similarities to our 
model. They investigated the properties of dark matter that self-annihilated 
(that is, interactions occurred in which the annihilation products were the 
dark matter particles). Their model affects the structure of halos by heating 
dense regions, and affects fluctuations in the early universe through the intro- 
duction of a non-trivial Jean's length for the dark matter. Our model differs in 
that we assume the dark matter annihilates to something relativistic that no 
longer contributes to the gravitational potential of a halo after annihilation. 
This should lead to a greater effect on the core properties of the halo, through 
adiabatic changes to the gravitational potential, than the model of |Carlson 



et al. (1992| ). The interacting dark matter model proposed by [Spergel fc Stein 



hardt (2000|) is constructed so as to avoid excessive suppression of the power 
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spectrum on small scales. The dark matter is presumed to be cold, and the 
interaction cross sections are small enough that the probability of interaction 
is small at the time of recombination. 

There are additional problems with an annihilating dark matter model. First, 
the amount of energy released in dark matter annihilation is tremendous; in 
our simulation of a single halo, the mass loss from annihilation (in whatever 
the decay products are) has an energy equivalent of 4.7 x 10 45 erg s _1 averaged 
over the lifetime of the universe. This is far too large to be consistent with 
observations of the gamma ray background. If the dark matter annihilates, 
the decay products cannot include photons. Second, as discussed at length by 
Kaplinghat et al. (2000|) , the particle physics underlying an annihilating dark 
matter model must be carefully contrived so as to avoid complete annihilation 
of the dark matter in the early universe. The suppression of annihilation must 
be in place for z > 10000, which is a rather low energy for any reasonable 
phase transitions in the dark matter sector. 



5 Conclusion 

Annihilating dark matter leads to halo density profiles consistent with those 
observed in dwarf spiral galaxies. Annihilation leads to a rapid reduction in the 
central density of a dark matter halo, driven mainly by adiabatic expansion 
of the core region as matter is lost. An annihilation cross section sufficient to 
generate a core in a dwarf galaxy halo generates a moderate core in a Milky 
Way sized halo, and a negligible core in a rich cluster. 

The core radius decreases rapidly with halo mass, primarily because larger 
objects form later in hierarchical clustering models, when the mean density 
of the universe is lower. Objects smaller than ~ 3 x 1O 8 M are unlikely to 
be present in an annihilating dark matter model, because the core radius will 
be large compared to the virial radius. A reduction in structure on this scale 
would be useful; CDM predicts many more halos of this size than are observed. 

On the other hand, there is no compelling motivation for such a complex dark 
matter candidate. A special mechanism is needed to suppress annihilation 
of the dark matter in the early universe, and the dark matter must decay 
to particles that do not interact with ordinary matter. Although a particle 
physics model with the required properties could no doubt be constructed, 
such a model would be contrived. 

Further study of the impact of annihilation in a cosmological setting would be 
useful in evaluating the promise of this model. 
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